## pval_cutoff: 0.05
## lfc_cutoff: 1
## low_counts_cutoff: 10

General statistics

# Number of samples
length(counts_data)
## [1] 6
# Number of genes
nrow(counts_data)
## [1] 55487
# Total counts
colSums(counts_data)
## SRR13535276 SRR13535278 SRR13535280 SRR13535300 SRR13535302 SRR13535304 
##     3107284     2321609     3701956     2491487     1580539     1861995

Create DDS objects

# Create DESeqDataSet object
dds <- get_DESeqDataSet_obj(counts_data, ~ treatment)
## [1] TRUE
## [1] TRUE
## [1] "DESeqDataSet object of length 55487 with 0 metadata columns"
## [1] "DESeqDataSet object of length 14648 with 0 metadata columns"
colData(dds)
## DataFrame with 6 rows and 25 columns
##              Assay Type AvgSpotLen       Bases  BioProject    BioSample      Bytes Center Name     Consent DATASTORE filetype DATASTORE provider               DATASTORE region  Experiment treatment GEO_Accession (exp)          Instrument LibraryLayout LibrarySelection  LibrarySource     Organism    Platform                    label ReleaseDate Sample Name                   source_name   SRA Study
##             <character>  <numeric>   <numeric> <character>  <character>  <numeric> <character> <character>        <character>        <character>                    <character> <character>  <factor>         <character>         <character>   <character>      <character>    <character>  <character> <character>                 <factor>   <POSIXct> <character>                   <character> <character>
## SRR13535276     RNA-Seq        300  8225466000 PRJNA694971 SAMN17588686 3252113587         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943360         A          GSM5043430 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043430 C2C12 proliferating myoblasts   SRP303354
## SRR13535278     RNA-Seq        300  9203426700 PRJNA694971 SAMN17588684 3619152333         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943362         A          GSM5043433 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043433 C2C12 proliferating myoblasts   SRP303354
## SRR13535280     RNA-Seq        300  9323939700 PRJNA694971 SAMN17588682 3735905901         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943364         A          GSM5043436 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043436 C2C12 proliferating myoblasts   SRP303354
## SRR13535300     RNA-Seq        300 12820015200 PRJNA694971 SAMN17587361 5047533646         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943384         E          GSM5043471 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA                  on land  2021-09-09  GSM5043471 C2C12 proliferating myoblasts   SRP303354
## SRR13535302     RNA-Seq        300 12499917600 PRJNA694971 SAMN17587359 4941074444         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943386         E          GSM5043475 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA                  on land  2021-09-09  GSM5043475 C2C12 proliferating myoblasts   SRP303354
## SRR13535304     RNA-Seq        300  7150086300 PRJNA694971 SAMN17587357 2845819297         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943388         E          GSM5043478 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA                  on land  2021-09-09  GSM5043478 C2C12 proliferating myoblasts   SRP303354

Sample-to-sample comparisons

# Transform data (blinded rlog)
rld <- get_transformed_data(dds)

PCA plot

pca <- rld$pca
pca_df <- cbind(as.data.frame(colData(dds)) %>% rownames_to_column(var = 'name'), pca$x)
summary(pca)
## Importance of components:
##                            PC1     PC2     PC3      PC4      PC5       PC6
## Standard deviation     40.9653 35.5987 26.4305 19.10170 17.10922 1.017e-13
## Proportion of Variance  0.3901  0.2946  0.1624  0.08482  0.06805 0.000e+00
## Cumulative Proportion   0.3901  0.6847  0.8471  0.93195  1.00000 1.000e+00
ggplot(pca_df, aes(x = PC1, y = PC2, color = label)) +
  geom_point() +
  geom_text(aes(label = name), position = position_nudge(y = -2), show.legend = F, size = 3) +
  scale_color_manual(values = colors_default) +
  scale_x_continuous(expand = c(0.2, 0))

Correlation heatmap

pheatmap(
  cor(rld$matrix),
  annotation_col = as.data.frame(colData(dds)) %>% select(label),
  color = brewer.pal(8, 'YlOrRd')
)

Wald test results

# DE analysis using Wald test
dds_full <- DESeq(dds)
colData(dds_full)
## DataFrame with 6 rows and 26 columns
##              Assay Type AvgSpotLen       Bases  BioProject    BioSample      Bytes Center Name     Consent DATASTORE filetype DATASTORE provider               DATASTORE region  Experiment treatment GEO_Accession (exp)          Instrument LibraryLayout LibrarySelection  LibrarySource     Organism    Platform                    label ReleaseDate Sample Name                   source_name   SRA Study        sizeFactor
##             <character>  <numeric>   <numeric> <character>  <character>  <numeric> <character> <character>        <character>        <character>                    <character> <character>  <factor>         <character>         <character>   <character>      <character>    <character>  <character> <character>                 <factor>   <POSIXct> <character>                   <character> <character>         <numeric>
## SRR13535276     RNA-Seq        300  8225466000 PRJNA694971 SAMN17588686 3252113587         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943360         A          GSM5043430 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043430 C2C12 proliferating myoblasts   SRP303354 0.970042620530961
## SRR13535278     RNA-Seq        300  9203426700 PRJNA694971 SAMN17588684 3619152333         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943362         A          GSM5043433 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043433 C2C12 proliferating myoblasts   SRP303354  1.33983479284256
## SRR13535280     RNA-Seq        300  9323939700 PRJNA694971 SAMN17588682 3735905901         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943364         A          GSM5043436 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043436 C2C12 proliferating myoblasts   SRP303354  1.08830437727923
## SRR13535300     RNA-Seq        300 12820015200 PRJNA694971 SAMN17587361 5047533646         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943384         E          GSM5043471 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA                  on land  2021-09-09  GSM5043471 C2C12 proliferating myoblasts   SRP303354  1.43563336306498
## SRR13535302     RNA-Seq        300 12499917600 PRJNA694971 SAMN17587359 4941074444         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943386         E          GSM5043475 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA                  on land  2021-09-09  GSM5043475 C2C12 proliferating myoblasts   SRP303354 0.769364687058128
## SRR13535304     RNA-Seq        300  7150086300 PRJNA694971 SAMN17587357 2845819297         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3.us-east-1  SRX9943388         E          GSM5043478 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA                  on land  2021-09-09  GSM5043478 C2C12 proliferating myoblasts   SRP303354 0.596304865135331
# Wald test results
res <- results(
  dds_full,
  contrast = c('treatment', condition, control),
  alpha = pval_cutoff
)
res
## log2 fold change (MLE): treatment A vs E 
## Wald test p-value: treatment A vs E 
## DataFrame with 14648 rows and 6 columns
##                            baseMean        log2FoldChange             lfcSE                 stat             pvalue              padj
##                           <numeric>             <numeric>         <numeric>            <numeric>          <numeric>         <numeric>
## ENSMUSG00000025900   4.914147077455     -4.92038939858524  2.01247166836637    -2.44494840644359 0.0144872864884619                NA
## ENSMUSG00000098104 4.09533781074173     0.853116865234774  1.10703456347057    0.770632546973277  0.440924763843958                NA
## ENSMUSG00000033845 107.622165011027    -0.107664037289825 0.423913828571251   -0.253976232982758  0.799513919532903 0.928640368582514
## ENSMUSG00000102275 2.36352235488834    -0.391339522802433  1.46485014416748   -0.267153281419679  0.789351145202682                NA
## ENSMUSG00000025903 97.3741809067814 -0.000670710765897874 0.485955321775427 -0.00138019018589496   0.99889876790933 0.999586418726237
## ...                             ...                   ...               ...                  ...                ...               ...
## ENSMUSG00000061654 1.69274896880504      1.42421218490796  2.57607432131434    0.552861450123578   0.58035828636321                NA
## ENSMUSG00000079834 28.8069677321529     0.998795513949071  0.80844927148458     1.23544611786829  0.216664518073648 0.549180345002018
## ENSMUSG00000095041 184.206277782681    0.0979395592239719 0.573810611348214    0.170682725775766  0.864473246335269 0.954515887544127
## ENSMUSG00000063897 31.5444997848201     -0.20997393171239 0.718744861367948   -0.292139732745717  0.770179788725231 0.915940742518338
## ENSMUSG00000095742 10.1107048409608      0.14999961155862 0.948601163242324    0.158127163839775  0.874356596236297 0.957516437938579
mcols(res)
## DataFrame with 6 rows and 2 columns
##                        type                               description
##                 <character>                               <character>
## baseMean       intermediate mean of normalized counts for all samples
## log2FoldChange      results  log2 fold change (MLE): treatment A vs E
## lfcSE               results          standard error: treatment A vs E
## stat                results          Wald statistic: treatment A vs E
## pvalue              results       Wald test p-value: treatment A vs E
## padj                results                      BH adjusted p-values
summary(res)
## 
## out of 14648 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 312, 2.1%
## LFC < 0 (down)     : 184, 1.3%
## outliers [1]       : 179, 1.2%
## low counts [2]     : 2840, 19%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plotDispEsts(dds_full)

Summary details

# Upregulated genes (LFC > 0)
res_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res[which(is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs E 
## Wald test p-value: treatment A vs E 
## DataFrame with 179 rows and 6 columns
##                            baseMean     log2FoldChange             lfcSE               stat    pvalue      padj
##                           <numeric>          <numeric>         <numeric>          <numeric> <numeric> <numeric>
## ENSMUSG00000067780  318.83564782775   -2.3706959444743  1.42869895255581  -1.65933903726418        NA        NA
## ENSMUSG00000025981 152.077989616933 -0.293486203907577  1.15037203790121  -0.25512285959508        NA        NA
## ENSMUSG00000038349 100.799409633029  -3.94260165031521  1.21130442878514   -3.2548396229917        NA        NA
## ENSMUSG00000026024 50.6697012550001  -3.51747879092895  1.25525049232079  -2.80221263600193        NA        NA
## ENSMUSG00000085842 21.9171266358604   5.16454652547104  2.11409207622401    2.4429146599402        NA        NA
## ...                             ...                ...               ...                ...       ...       ...
## ENSMUSG00000005871 403.634302628017 -0.524237280126561 0.979145628125162 -0.535402768565034        NA        NA
## ENSMUSG00000044595 41.1231231533217   1.83318093904902   1.6810418210932   1.09050287509021        NA        NA
## ENSMUSG00000024597  346.21374732627  -1.29473000804499 0.973881588495617  -1.32945321416846        NA        NA
## ENSMUSG00000118138  23.291043050614   5.90635654386202   3.5290227276766   1.67365216935018        NA        NA
## ENSMUSG00000033417 290.918339858813 -0.928829346925173 0.966631349728102 -0.960893051095888        NA        NA
# Low counts (only padj is NA)
res[which(is.na(res$padj) & !is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs E 
## Wald test p-value: treatment A vs E 
## DataFrame with 2840 rows and 6 columns
##                            baseMean     log2FoldChange            lfcSE               stat             pvalue      padj
##                           <numeric>          <numeric>        <numeric>          <numeric>          <numeric> <numeric>
## ENSMUSG00000025900   4.914147077455  -4.92038939858524 2.01247166836637  -2.44494840644359 0.0144872864884619        NA
## ENSMUSG00000098104 4.09533781074173  0.853116865234774 1.10703456347057  0.770632546973277  0.440924763843958        NA
## ENSMUSG00000102275 2.36352235488834 -0.391339522802433 1.46485014416748 -0.267153281419679  0.789351145202682        NA
## ENSMUSG00000102135  4.9457418399281 -0.328110499112732 1.08089829259967 -0.303553536312462  0.761468054712299        NA
## ENSMUSG00000098201 2.01346784837404 -0.896309374862552 1.72861159166932 -0.518514037035346  0.604099668891218        NA
## ...                             ...                ...              ...                ...                ...       ...
## ENSMUSG00000064344 2.73922713299793  0.163917921084994 1.41591685542737  0.115768041362512  0.907836379470033        NA
## ENSMUSG00000064349 3.00782467550163  -0.12734120960026 1.24974955476499 -0.101893382649939  0.918841302505899        NA
## ENSMUSG00000064358 2.70134753598084  0.141492229313294 1.65446199795407 0.0855215952305128   0.93184672785224        NA
## ENSMUSG00000064369 4.23154180234235   1.42531627839752 1.24382109968624   1.14591743037409  0.251829318362109        NA
## ENSMUSG00000061654 1.69274896880504   1.42421218490796 2.57607432131434  0.552861450123578   0.58035828636321        NA

Shrunken LFC results

plotMA(res)

# Shrunken LFC results
res_shrunken <- lfcShrink(
  dds_full,
  coef = str_c('treatment_', condition, '_vs_', control),
  type = 'apeglm'
)
res_shrunken
## log2 fold change (MAP): treatment A vs E 
## Wald test p-value: treatment A vs E 
## DataFrame with 14648 rows and 5 columns
##                            baseMean       log2FoldChange             lfcSE             pvalue              padj
##                           <numeric>            <numeric>         <numeric>          <numeric>         <numeric>
## ENSMUSG00000025900   4.914147077455   -0.241445095152026 0.566057167536038 0.0144872864884619 0.152194512661901
## ENSMUSG00000098104 4.09533781074173    0.140171540467965 0.461172330088417  0.440924763843958 0.746711563326519
## ENSMUSG00000033845 107.622165011027  -0.0595416997344277 0.318135895241512  0.799513919532903  0.92952733548211
## ENSMUSG00000102275 2.36352235488834  -0.0364245117406618 0.454994221479306  0.789351145202682                NA
## ENSMUSG00000025903 97.3741809067814 0.000553662408882452 0.339890271773017   0.99889876790933 0.999554374615645
## ...                             ...                  ...               ...                ...               ...
## ENSMUSG00000061654 1.69274896880504   0.0466582594287195 0.472044184129129   0.58035828636321                NA
## ENSMUSG00000079834 28.8069677321529    0.285902217369414 0.497015110137123  0.216664518073648  0.55206916045031
## ENSMUSG00000095041 184.206277782681   0.0410310477560653 0.367468135932862  0.864473246335269 0.954844930608897
## ENSMUSG00000063897 31.5444997848201  -0.0639566027984351 0.400050017030481  0.770179788725231 0.916738662708078
## ENSMUSG00000095742 10.1107048409608   0.0313236418014735 0.425676345462796  0.874356596236297 0.958102496428203
plotMA(res_shrunken)

mcols(res_shrunken)
## DataFrame with 5 rows and 2 columns
##                        type                               description
##                 <character>                               <character>
## baseMean       intermediate mean of normalized counts for all samples
## log2FoldChange      results  log2 fold change (MAP): treatment A vs E
## lfcSE               results            posterior SD: treatment A vs E
## pvalue              results       Wald test p-value: treatment A vs E
## padj                results                      BH adjusted p-values
summary(res_shrunken, alpha = pval_cutoff)
## 
## out of 14648 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 308, 2.1%
## LFC < 0 (down)     : 175, 1.2%
## outliers [1]       : 179, 1.2%
## low counts [2]     : 2272, 16%
## (mean count < 4)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Summary details

# Upregulated genes (LFC > 0)
res_shrunken_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_shrunken_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res_shrunken[which(is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs E 
## Wald test p-value: treatment A vs E 
## DataFrame with 179 rows and 5 columns
##                            baseMean      log2FoldChange             lfcSE    pvalue      padj
##                           <numeric>           <numeric>         <numeric> <numeric> <numeric>
## ENSMUSG00000067780  318.83564782775   -0.22411703750797 0.536162264989718        NA        NA
## ENSMUSG00000025981 152.077989616933 -0.0416595659544155 0.442086642338937        NA        NA
## ENSMUSG00000038349 100.799409633029   -2.96028380309171  1.53835977480067        NA        NA
## ENSMUSG00000026024 50.6697012550001   -0.55794504867263  1.14073581036123        NA        NA
## ENSMUSG00000085842 21.9171266358604   0.171735351785563 0.521147556380492        NA        NA
## ...                             ...                 ...               ...       ...       ...
## ENSMUSG00000005871 403.634302628017 -0.0999653241456368 0.440350464820972        NA        NA
## ENSMUSG00000044595 41.1231231533217   0.129395103559718 0.487081495732831        NA        NA
## ENSMUSG00000024597  346.21374732627  -0.271004345889595 0.521619060100433        NA        NA
## ENSMUSG00000118138  23.291043050614  0.0577703753824257 0.480736973051452        NA        NA
## ENSMUSG00000033417 290.918339858813  -0.187614640240523 0.470764277170404        NA        NA
# Low counts (only padj is NA)
res_shrunken[which(is.na(res_shrunken$padj) & !is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs E 
## Wald test p-value: treatment A vs E 
## DataFrame with 2272 rows and 5 columns
##                            baseMean      log2FoldChange             lfcSE            pvalue      padj
##                           <numeric>           <numeric>         <numeric>         <numeric> <numeric>
## ENSMUSG00000102275 2.36352235488834 -0.0364245117406618 0.454994221479306 0.789351145202682        NA
## ENSMUSG00000098201 2.01346784837404 -0.0631972789924058 0.464941327316382 0.604099668891218        NA
## ENSMUSG00000103903 3.35339380940403  0.0953761002836069 0.477444385207544 0.365967184382433        NA
## ENSMUSG00000079671 1.77247301755752  -0.031554464222372 0.459385323682869 0.799750938797838        NA
## ENSMUSG00000083422 2.09976591294719  -0.118324077380156 0.484283472973611  0.29066750440083        NA
## ...                             ...                 ...               ...               ...       ...
## ENSMUSG00000064342 2.42847689971841  0.0625177992571637 0.469762522277571 0.546404966544711        NA
## ENSMUSG00000064344 2.73922713299793  0.0168427484873255 0.451689823183754 0.907836379470033        NA
## ENSMUSG00000064349 3.00782467550163 -0.0158356528673377 0.444956540472799 0.918841302505899        NA
## ENSMUSG00000064358 2.70134753598084   0.010794852575915 0.458128601167714  0.93184672785224        NA
## ENSMUSG00000061654 1.69274896880504  0.0466582594287195 0.472044184129129  0.58035828636321        NA

Visualizing results

Heatmaps

# Plot normalized counts (z-scores)
pheatmap(counts_sig_norm[2:7], 
         color = brewer.pal(8, 'YlOrRd'), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label),
         border_color = NA,
         fontsize = 10,
         scale = 'row',
         fontsize_row = 10, 
         height = 20)

# Plot log-transformed counts
pheatmap(counts_sig_log[2:7], 
         color = rev(brewer.pal(8, 'RdYlBu')), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label),
         border_color = NA,
         fontsize = 10,
         fontsize_row = 10, 
         height = 20)

# Plot log-transformed counts (top 24 DE genes)
pheatmap((counts_sig_log %>% filter(ensembl_gene_id %in% res_sig_df$ensembl_gene_id[1:24]))[2:7],
         color = rev(brewer.pal(8, 'RdYlBu')), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label), 
         fontsize = 10,
         fontsize_row = 10, 
         height = 20)

Volcano plots

# Unshrunken LFC
res_df %>% 
  mutate(
    sig_threshold = if_else(
      padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
      if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
      'non-DE'
    )
  ) %>% 
  filter(!is.na(sig_threshold)) %>% 
  ggplot() +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
  scale_color_manual(values = c('blue', 'red', 'gray')) +
  xlab('log2 fold change') + 
  ylab('-log10 adjusted p-value')

# Shrunken LFC
res_shrunken_df %>% 
  mutate(
    sig_threshold = if_else(
      padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
      if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
      'non-DE'
    )
  ) %>% 
  filter(!is.na(sig_threshold)) %>% 
  ggplot() +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
  scale_color_manual(values = c('blue', 'red', 'gray')) +
  xlab('log2 fold change') + 
  ylab('-log10 adjusted p-value')

GSEA (all)

Hallmark genesets

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

GSEA (DE)

Hallmark genesets

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

System Info

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] fgsea_1.12.0                Rcpp_1.0.3                  RColorBrewer_1.1-2          pheatmap_1.0.12             DESeq2_1.26.0               SummarizedExperiment_1.16.1 DelayedArray_0.12.3         BiocParallel_1.20.1         matrixStats_0.57.0          Biobase_2.46.0              GenomicRanges_1.38.0        GenomeInfoDb_1.22.1         IRanges_2.20.2              S4Vectors_0.24.4            BiocGenerics_0.32.0         scales_1.1.1                forcats_0.4.0               stringr_1.4.0               dplyr_1.0.2                 purrr_0.3.3                 readr_1.3.1                 tidyr_1.0.0                 tibble_3.1.0                ggplot2_3.3.3               tidyverse_1.2.1            
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1       ellipsis_0.3.0         htmlTable_1.13.3       XVector_0.26.0         base64enc_0.1-3        rstudioapi_0.10        farver_2.1.0           bit64_0.9-7            mvtnorm_1.1-1          apeglm_1.8.0           AnnotationDbi_1.48.0   fansi_0.4.0            lubridate_1.7.4        xml2_1.2.2             splines_3.6.3          geneplotter_1.64.0     knitr_1.25             Formula_1.2-3          jsonlite_1.6           broom_0.7.5            annotate_1.64.0        cluster_2.1.0          png_0.1-7              compiler_3.6.3         httr_1.4.1             backports_1.1.5        assertthat_0.2.1       Matrix_1.2-18          cli_1.1.0              acepack_1.4.1          htmltools_0.5.1.1      tools_3.6.3            coda_0.19-3            gtable_0.3.0           glue_1.4.2             GenomeInfoDbData_1.2.2 fastmatch_1.1-0        bbmle_1.0.23.1         cellranger_1.1.0       jquerylib_0.1.3        vctrs_0.3.4            xfun_0.22              rvest_0.3.5            lifecycle_0.2.0        XML_3.99-0.3           MASS_7.3-51.5          zlibbioc_1.32.0        hms_0.5.2              yaml_2.2.0             memoise_1.1.0          gridExtra_2.3          emdbook_1.3.12         sass_0.3.1             bdsmatrix_1.3-4        rpart_4.1-15           latticeExtra_0.6-29    stringi_1.4.3          RSQLite_2.2.1          genefilter_1.68.0      checkmate_1.9.4        rlang_0.4.8            pkgconfig_2.0.3        bitops_1.0-6           evaluate_0.14          lattice_0.20-38        labeling_0.3           htmlwidgets_1.5.1      bit_1.1-15.1           tidyselect_1.1.0       plyr_1.8.4             magrittr_1.5           R6_2.4.0               generics_0.0.2         Hmisc_4.3-0            DBI_1.1.0              pillar_1.5.1           haven_2.2.0            foreign_0.8-75         withr_2.1.2            survival_3.1-8         RCurl_1.95-4.12        nnet_7.3-12            modelr_0.1.5           crayon_1.3.4           utf8_1.1.4             rmarkdown_2.7          jpeg_0.1-8.1           locfit_1.5-9.4         grid_3.6.3             readxl_1.3.1           data.table_1.13.6      blob_1.2.1             digest_0.6.27          xtable_1.8-4           numDeriv_2016.8-1.1    munsell_0.5.0          bslib_0.2.4